# Replication archive for:
# Galos and Coppock. 2023. 
# Gender Composition Predicts Gender Bias: A Meta-reanalysis of Hiring Discrimination Audit Experiments. 
# Science Advances. Forthcoming.

library(tidyverse)
library(xtable)
library(metafor)
library(broom)
source("code/helpers.R")

occupation_cates <-
  read_rds("data/occupation_cates.rds")

study_ids <- unique(occupation_cates$study_id)

gender_gradients_list <- list()

for (i in seq_along(study_ids)) {
  local_id <- study_ids[i]
  
  local_dat <-
    occupation_cates |>
    filter(study_id == local_id) |>
    mutate(
      occupation_label = paste0(
        occupation_observed,
        "\n",
        round(fraction_female * 100, 0)
      ),
      entry = make_se_entry(estimate, std.error),
      facet = "estimates"
    )
  
  if (nrow(local_dat) > 2) {
    local_meta <-
      rma.uni(
        estimate,
        sei = std.error,
        mods = ~ fraction_female,
        data = local_dat
      ) |>
      tidy(conf.int = TRUE) |>
      filter(term == "fraction_female") |>
      mutate(
        occupation_label = factor("Gender gradient"),
        entry = make_se_entry(estimate, std.error),
        facet = "meta"
      )
    
    gender_gradients_list[[i]] <- local_meta |> 
      mutate(study_id = local_id)
    
    local_dat <-
      local_dat |>
      bind_rows(local_meta)
    
  }
  
  local_dat <-
    local_dat |>
    ungroup() |>
    mutate(occupation_label = fct_reorder(occupation_label, fraction_female))
  
  g1 <-
    ggplot(local_dat, aes(estimate, occupation_label)) +
    geom_vline(
      xintercept = 0,
      color = "red",
      linetype = "dashed",
      alpha = 0.5
    ) +
    geom_point() +
    geom_linerange(aes(xmin = conf.low, xmax = conf.high)) +
    geom_text(aes(label = entry), nudge_y = 0.3, size = 3) +
    facet_grid(rows = vars(facet),
               space = "free",
               scales = "free_y") +
    coord_cartesian(xlim = c(-100, 100)) +
    theme_bw() +
    theme(
      strip.text = element_blank(),
      text = element_text(size = 8),
      axis.title.y = element_blank()
    ) +
    labs(x = "CATE estimate (top facet) or gender gradient estimate (bottom facet)")
  
  # ggsave(
  #   paste0("for_overleaf/figures/", study_ids[i], "_occupation_cates.pdf"),
  #   plot = g1,
  #   width = 6.5,
  #   height = 1.5 + nrow(local_dat) * 0.25
  # )
}

# this outputs gender gradients for further analysis
gender_gradients_df <-
  gender_gradients_list |>
  bind_rows() 

write_rds(gender_gradients_df, "data/gender_gradients_df.rds")
